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We investigate the dynamics of stellar-mass black holes (BH) in star clusters focusing on the dy- 
namical formation of BH-BH binaries, which are very important sources of gravitational waves 
(GW). We examine the properties of these BH-BH binaries through direct N-body computations 
of Plummer clusters, having initially N(0) < 10 s low mass stars and a population of stellar mass 
BHs, using the state-of-the-art N-body integrator "NBODY6". We find that the stellar mass BHs 
segregate rapidly into the cluster core to form a central dense sub-cluster of BHs in which BH- 
BH binaries form via 3-body encounters. While most of the BH binaries finally escape from the 
cluster by recoils due to super-elastic encounters with the single BHs, we find that for clusters 
with N(0) ^ 5 x 10 4 , typically a few of them dynamically harden to the extent that they can merge 
via GW emission within the cluster. Also for each of such clusters, there are a few escaped BH 
binaries that merge within a Hubble time, most of the merger happening within a few Gyr of 
cluster evolution. These results imply that the intermediate-aged massive clusters constitute the 
most important class of star clusters that can produce dynamical BH-BH mergers at the present 
epoch. The BH-BH merger rates obtained from our computations imply a significant detection 
rate (w 30 yr~ 1 ) for the "Advanced LIGO" GW detector that will become operative in the near 
future. Finally, we briefly discuss our ongoing development on this work incorporating the for- 
mation of BHs in star clusters from stellar evolution. In particular, we highlight on the effect of 
stellar metallicity on the BH sub-cluster driven expansion of a star cluster's core. 
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1. Introduction 

Star clusters, e.g., globular clusters (henceforth GC), young and intermediate-age massive 
clusters harbor a large overdensity of compact stellar remnants compared to that in the field by 
virtue of their higher stellar density and the mass segregation of the stellar remnants. These com- 
pact stars, which are neutron stars (hereafter NS) and black holes (hereafter BH), are produced 
through the evolution of the most massive stars in the cluster within the first ~ 50 Myr after cluster 
formation. Since the BHs are significantly heavier than the remaining low-mass stars of the cluster, 
they undergo a runaway mass-segregation or mass-stratification instability [1] provided a signifi- 
cant fraction of them are retained in the cluster following their progenitor supernovae. Due to this 
instability (also known as the Spitzer instability), the BHs concentrate quickly to the cluster core 
forming a dense sub-cluster made purely of BHs which is nearly dynamically isolated from rest of 
the cluster. 

Earlier studies of the dynamical properties of such BH sub-clusters (e.g., references [2-4]) 
indicate that the density of the BHs there is large enough that the BH-BH binary formation through 
3-body encounters [5] becomes important. These dynamically formed BH binaries then shrink 
in orbit or "harden" through repeated super-elastic encounters with the surrounding BHs (Heggie's 
law [6]). The binding energy of the BH binaries released in this way is earned by the single BHs and 
BH binaries involved in the encounters. This causes the single and the binary BHs to get ejected 
from the BH core to the outer regions of the cluster and they deposit their energy to the cluster 
stars while returning to the core via dynamical friction [1, 5], thereby heating the cluster. This 
heating mechanism is most efficient in the central region of the cluster since the stellar density is 
there the highest. As the BH binaries harden, the encounter-driven recoils in the BH core become 
increasingly stronger and finally the recoils are large enough that the encountering single BHs 
and/or the BH binaries escape from the cluster. This also results in heating of the cluster core 
as the associated mass-loss dilutes the potential well there. These heating mechanisms result in 
an expansion of the central part of the cluster and modify the dynamical evolution of the cluster 
significantly. 

In the present work, we make a detailed study of the dynamics of the BH-BH binaries formed 
in a BH sub-cluster. Specifically, we investigate whether hard enough BH binaries, that can merge 
via gravitational radiation in a Hubble time within the cluster or after being ejected, can be formed 
in such a sub-cluster. To that end, we perform a series of direct N-body integrations of concentrated 
star clusters (half-mass radius r h < 1 pc) consisting of N ~ 10 4 — 10 5 low-mass stars, in which a 
certain number of stellar-mass BHs are added, representing a star cluster with an evolved stellar 
population. However, we also consider initial clusters with a full stellar mass spectrum in which 
the BHs are formed from stellar evolution. 

2. Computations 

To study the dynamics of BHs in star clusters, we perform direct N-body integrations of model 
star clusters using the "NBODY6" N-body integration code running on Graphical Processing Unit 
(GPU) platforms [7]. The initial clusters follow a Plummer model with half-mass radius 77,(0) < 1 
pc consisting of N(0) < 10 5 low-mass main-sequence stars in the mass-range O.5M < m < 1 .0M & . 
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Figure 1: Mass segregation of BHs shown by the radial position R vs. time t (left panel) for the model 
C50K80. The BHs segregate in ~ 50 Myr during which Nbh (the number of BHs bound to the cluster) 
remains constant (right panel). As the BH sub-cluster becomes dense enough that BH-BH binaries begin to 
form via 3-body encounters, single BHs and BH binaries start escaping from the cluster depleting Nbh- 



A specified number of BHs are added to each cluster with the same initial Plummer distribution as 
the stars. The initial number of BHs are chosen according to the BH retention fraction, where we 
explore both the full retention and the case where half of the BHs have escaped from the cluster by 
supernova natal kicks. For simplicity, we consider only equal-mass BHs in the present work, with 
the representative value of Mbh = 1QM©. Such an initial condition mimics the epoch at which the 
massive stars have already been evolved to produce their remnant BHs. 

To evolve the BH-BH binaries due to GW radiation, the well known Peters' formula is adapted 
in NBODY6, according to which the merger time T mrg of an equal-mass BH-BH binary due to 
gravitational wave (GW) emission is given by, 

T ^™M0(£^ i - e2)V2 - <21) 

where, a is the semi-major-axis of the binary and e is its eccentricity. Studies in numerical rel- 
ativity indicate that for unequal-mass BH-BH mergers or even for equal-mass BHs with unequal 
spins, the merged BH product acquires a velocity -kick of typically 100 km s _1 or more due to 
the asymmetric momentum outflow from the system, associated with the asymmetry of the GW 
emission. Therefore, in our computations, we provide an arbitrarily large velocity kick of 150 km 
s~ ! to the merged BH immediately after a BH-BH merger occurs within the cluster, to make sure 
that it escapes. Also, we evolve the clusters isolated, i.e., in absence of a galactic tidal field, as the 
formation of the BH-core through mass segregation and its dynamics remain largely unaffected by 
the presence of a tidal field, which mainly influences the stars near the tidal boundary. Furthermore, 
for simplicity, we do not consider primordial binaries in this initial study. See reference [4] for the 
details on the implementation of initial conditions and GW natal kicks. 

2.1 Mimicking a GC core: reflective boundary 

We also perform computations with a smaller number of stars and BHs that are confined within 
a reflecting spherical boundary [4]. With such a dynamical system, one can mimic the core of a 
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Table 1: Summary of the computations for isolated clusters and those with reflective boundary. The different 
columns represent as follows: Col. (1): identity of the particular model — similar values with different 
names (ending with A,B etc) imply runs repeated with different random seeds, Col. (2): total number of 
stars N(Q) initially, Col. (3): number of computations N s j,„ with the particular cluster, Col. (4): initial half- 
mass radius of the cluster 0,(0) (isolated cluster) or radius of reflective sphere R s , Col. (5): initial number of 
BHs Nbh(0), Col. (6): total number of BH-BH binary mergers within the cluster N mrg , Col. (7): the times 
t mrg corresponding to the mergers, Col. (8): number of escaped BH-pairs N esc — the three values of N esc are 
those with T mrg < 3 Gyr, 1 Gyr and 100 Myr respectively. 



Model name 


2V(0) 




r h (0) oiR s (pc) 


N BH (P) 


^Yi mrg 


t mrg (Myr) 


Nesc 


Isolated clusters 


C10K20 


10000 


10 


1.0 


20 









C25K50 


25000 


10 


1.0 


50 







3 1 1 


C50K80 


45000 




1.0 


80 


1 


698.3 


3 1 


C50K80.1 


45000 




0.5 


80 


2 


217.1,236.6 


32 1 


C50K40.1 


45000 




0.5 


40 







1 1 1 


C50K200 


50000 




1.0 


200 


2 


100.8, 467.8 


000 


C65K110 


65000 




1.0 


110 


1 


314.6 


42 1 


C65K110.1 


65000 




0.5 


110 







43 1 


C65K55.1 


65000 




0.5 


55 


1 


160.5 


100 


C100K80 


100000 




1.0 


80 


2 


219.4, 603.2 


521 


C100K200 


100000 




1.0 


200 







544 


Reflective boundary 


R3K180 


3000 




0.4 


180 


1 


1723.9 


5 3 1 


R4K180A 


4000 




0.4 


180 


1 


3008.8 


22 1 


R4K180B 


4000 




0.4 


180 


2 


100.2, 1966.5 


2 1 


R3K100 


3000 




0.4 


100 


2 


3052.8, 3645.9 


1 1 


R4K100A 


4000 




0.4 


100 


2 


104.4, 814.2 


3 3 1 


R4K100B 


4000 




0.4 


100 


1 


1135.3 


3 3 3 



massive cluster, where the BHs are concentrated after mass segregation. The advantage of this 
approach is that one can compute the evolution of a massive cluster with many fewer stars, allowing 
much faster computation. We integrate N = 3000 — 4000 stars bounded within 0.4 pc which gives 
a stellar density of ~ 1O 4 M pc~ 3 appropriate for the core-density of a massive cluster. However, 
stars and BHs faster than a pre-assigned speed v esc ~ 24 km sec -1 , representing the escape speed 
from the parent cluster, are allowed to escape through the reflective boundary. 

3. Results 

Table 1 summarizes the results of our computations. For demonstration purposes, we choose 
one of the models, viz. , C50K80 (see Table 1 ) — other models generally possess similar properties. 
Fig. 1 (left panel) demonstrates the mass segregation of the BHs in the cluster which takes about 
50 Myr. As the BHs segregate within about 0.3 pc of the cluster core, the BH density of becomes 
large enough to form BH-BH binaries through 3-body encounters. Once BH-BH binaries begin 
forming, single BHs and BH binaries start escaping from the BH-core due to the encounter-recoils 
(see Sec. 1). In Fig. 1 (left panel), one can clearly distinguish the two phases of the BH subsystem 
— the initial segregation phase and the BH-core formation, the radial positions of the BHs being 
scattered outwards in the latter phase due to the recoils. The corresponding decrease of Nbh during 
this phase is also shown in Fig. 1 (right panel). 
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Figure 2: The positions of all the BH binaries within the cluster on the 1 — e 2 vs. a plane where the different 
symbols are used to distinguish between different BH pairs (for model C50K80). The color-coding of the 
points, as indicated by the color-scale, represents the time (in Myr) at which they appear at a particular 
location in the above plane. The lines of constant GW merger times, T mrg = constant, are overplotted. 

3.1 Mergers and escapers 

To study the possibility of BH-BH mergers, we consider the positions of the BH-BH binaries 
within the cluster in an a vs. (1 — e 2 ) plane as shown in Fig. 2, where each pair is represented 
by a different symbol and are color-coded with the cluster evolution time in Myr. For each BH- 
pair, a and e fluctuate over the plane, the changes occurring over a collision time-scale which is 
~ Myr. Since the orbital periods of the corresponding binaries are much shorter (from ~ days to 
years), these points generally represent binaries which are stable over many orbits. Overplotted in 
Fig. 2 are the lines of constant GW merger time T mrg as given by Eqn. (2.1). While most of the 
points correspond to very long merger times, a few of them do lie close to the T mrg = 10 Myr line. 
This implies that these binaries are indeed hardened up to small enough a and/or acquire sufficient 
eccentricity that had they been left unperturbed, they would have merged via GW emission within 
several Myr. However, these merging BH-pairs can still be perturbed by further encounters on 
timescales ~ Myr which can often prevent them from merging. In the particular example shown 
in Fig. 2, only one of them could merge within the cluster. This feature is found to be generally 
true for all the models reported here (see Table 1), i.e., each of them does produce BH-pairs that 
are capable of merging within several Myr. On the other hand, among the escaped BH binaries, 
all of them with GW merger times smaller than a Hubble time are of interest since they remain 
unperturbed afterwards. For each of our computed models, a few of the escaped BH-BH binaries 
can merge in T mrg < 3 Gyr as can be read from Table 1 (also, see reference [4]). We also find that 
for all of the computed models, most mergers occur within a few Gyr (see reference [4]). 

The above results imply that an important class of star clusters for dynamically forming BH 
binaries that merge at the present epoch are those with initial mass M c/ (0)>3 x 1O 4 M and are 
few Gyr old. Such clusters represent intermediate-age massive clusters (hereafter IMC) with initial 
masses close to the upper-limit of the initial cluster mass function (ICMF) in spiral and starburst 
galaxies [8]. Globular clusters, which are at least an order of magnitude more massive, are typi- 
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cally much older (~ 10 Gyr), so that they cannot contribute significantly to the present-day merger 
rate, as most of the mergers from them would have occurred earlier. On the other hand, young 
massive clusters, with ages < 50 Myr, are generally too young to produce BH-BH mergers, as the 
stratification of the BHs and the formation of the BH core itself takes longer. Hence, IMCs seem 
to be the most likely star clusters for dynamically generating present-day BH-BH merger events. 

3.2 Detection rate of BH-BH mergers 

The results from Table 1 give an average merger rate of ps 0.4/cluster/Gyr for those occurring 
within the clusters and ps 0.9/cluster/Gyr for the escapers, totaling S% mrg ~ 1.3/cluster/Gyr (see 
reference [4]). To estimate the detection rate of BH-BH mergers from IMCs, we adopt the space 
density, as derived by reference [9], for young populous clusters which is p c / ps 3.5 h 3 Mpc~ 3 , 
where h = 0.73 is the Hubble parameter. This merger rate implies that the IMCs would yield a BH- 
BH merger detection rate of ^AdLioo ~ 30 yr _1 for the upcoming "Advanced LIGO" (AdLIGO) 
GW detector (see reference [4] and citations therein). Interestingly, the above dynamical BH-BH 
merger detection rate can be more than 10 times higher than that from individual primordial binaries 
[10]. Our results then imply that the dynamically formed BH-BH binaries constitute the dominant 
contribution to the BH-BH merger detection from the Universe [4]. 

4. Introduction of stellar evolution 

Although the study described above has provided significant insights into the dynamics of 
stellar mass BHs in star clusters and its role in dynamically generating BH-BH merger events and 
has been performed self-consistently using direct N-body computations, the cluster models would 
be more realistic if the BHs are formed with a mass function from stellar evolution. To materialize 
such improvements, we model the initial Plummer clusters with all stars at their Zero- Age-Main- 
Sequence (ZAMS) having a Kroupa IMF [1 1] (upto IOOMq). These models are also evolved using 
the NBODY6 integrator where the direct N-body integration of the cluster stars is coupled with 
their stellar evolution using analytic but well tested stellar evolution recipes [12]. 

Fig. 3 shows the evolution of the core radii r c for two cluster models with different metallici- 
ties, viz. , Z = Z and Z = 0.01Z Q , having otherwise exactly the same initial conditions (N(0) = 10 5 , 
r/j (0) = 1 pc Plummer spheres and all stars having identical positions, velocities and masses fol- 
lowing a Kroupa IMF). All the BHs are retained in the clusters following their formation via super- 
novae. For the cluster with lower Z, the BH masses are considerably larger than that for the higher 
Z case (average BH mass p=i 14M for Z = 0.01Z Q as opposed to « 7M Q for Z = Z ) so that the ex- 
pansion of the core radius due to the heating effect of the BH sub-cluster (see Sec. 1) is significantly 
more prominent for the lower Z model, as can be seen in Fig. 3. Also, the Z = 0.01Z Q model have 
produced 3 BH-BH binaries among its escapers that can merge via GW radiation within a Hubble 
time while its high Z counterpart have produced only one. This can be expected as the heavier BHs 
in the lower Z cluster create a denser BH core where the binary BH formation rate through 3-body 
encounters and their subsequent dynamical hardening rates are higher. This indicates that the mod- 
ification of the dynamical evolution of a cluster by its BHs and the dynamical BH-BH merger rate 
depend crucially on its metallicity. Studying the effect of the stellar content of a cluster and of its 
properties on the dynamics of its BHs is under progress. 
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Figure 3: Evolution of core radii r c for two computed clusters with metallicities Z = Z Q and Z = 0.01Z E , 
but otherwise identical initial conditions (see text). A full BH retention is assumed. Such comparison 
demonstrates that the heating effect of the BH core and the resulting expansion of r c is significantly larger 
for the lower metallicity cluster (left panel) as the BHs formed are considerably more massive in that case. 
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